R Notebook

This notebook contains code for a homework problem for Stat 242 : Time Series in Fall 2024.

Original data come from

Fisher, N.I. (1993) Statistical Analysis of Circular Data. Cambridge: Cambridge University Press.

In the class, we modeled these data as coming from a Hidden Markov Model with \(Y_t\) measured wind direction viewed as noisy measurements of a true underlying wind direction \(\alpha_t\). The processed data in this repo are the result of a particle filter that I implement on these data (I omit this in case this is a homework question in future years of the course). After obtaining posterior samples for \(\alpha_t\) and using \(\mod(\alpha_t,2\pi)\) to restrict them to \((0,2\pi]\), I calculate the circular mean \(\beta_t\) on them (source) and then plot (\(\cos(\beta_t),\sin(\beta_t)\)) in the animation.

Version 1: This just plots the circular mean, 20 at a time. The 20 at a time is a bit strange-looking, but it does make it easier to visually track the trend in the motion over time than if I just showed a single dot. I “pad” each segment of 20 so that you see times 1-20, 20-40, 40-60 etc. with overlap (helps continuity of the motion).

Tmax = dim(df)[1]
t = 20
n_frames = Tmax/t 

# make each transition point appear double 
# repeat every t^th twice so e.g. if t=20, get 1:20, 20:40 etc.
ind = sapply(1:n_frames, function(i){(t*(i-1)):(t*(i-1)+t)})
ind = ind[2:length(ind)] #remove extra 0
dfpad = df[ind,]

# create frames
dfpad$frame =  c(rep(1,t),rep(2:n_frames, each = t+1))

plotly::plot_ly(data = dfpad,
                 x = ~eastwest, #~ so looks for them in dataset
                 y = ~northsouth,
                 frame = ~frame,
                 type = "scatter",
                 text = ~paste("Time:", Time),
                 mode = "lines+markers",
                 marker = list(size = 8,
                               symbol = "circle",
                               sizemode = "diameter"),
                 line = list(shape = "linear", width = 2)
                ) %>%
plotly::layout(xaxis = list(title = "East-West Direction"),
                 yaxis = list(title = "North-South Direction"),
                showlegend = F
        ) %>%
plotly::animation_opts(frame = 200,
                       transition = 10,
                       redraw = F) 

Version 2: This plots the circular mean and the 5% and 95% quantile of the posterior circular means, 5 at a time. It’s a bit funky to look at but it does give a better picture because it shows that the uncertainty in the wind direction is often quite large.

Tmax = dim(df)[1]
t = 5
n_frames = Tmax/t 

# add the lower and upper 
main = cbind(df[,c("Time","Wind","eastwest", "northsouth")], type = rep("mean",Tmax))
lower = cbind(df[,c("Time","Wind","eastwest.lower", "northsouth.lower")], type = rep("lower", Tmax))
upper = cbind(df[,c("Time","Wind","eastwest.upper", "northsouth.upper")], type = rep("upper", Tmax))
lower = rename(lower, "northsouth" = "northsouth.lower")
lower = rename(lower, "eastwest" = "eastwest.lower")
upper = rename(upper, "northsouth" = "northsouth.upper")
upper = rename(upper, "eastwest" = "eastwest.upper")


# Get pad indices
# make each transition point appear double 
# repeat every t^th twice so e.g. if t=20, get 1:20, 20:40 etc.
ind = sapply(1:n_frames, function(i){(t*(i-1)):(t*(i-1)+t)})
ind = ind[2:length(ind)] #remove extra 0
main = main[ind,]
lower = lower[ind,]
upper = upper[ind,]

# create frames and add to each
frame =  c(rep(1,t),rep(2:n_frames, each = t+1))
main$frame = frame
lower$frame = frame
upper$frame = frame

#consolidate
df_all = rbind(main, lower, upper)


plotly::plot_ly(data = df_all,
                 x = ~eastwest, #~ so looks for them in dataset
                 y = ~northsouth,
                 frame = ~frame,
                 color = ~type,
                 type = "scatter",
                 text = ~paste("Time:", Time),
                 mode = "lines+markers",
                 marker = list(size = 8,
                               symbol = "circle",
                               sizemode = "diameter"),
                 line = list(shape = "linear", width = 2)
                ) %>%
plotly::layout(xaxis = list(title = "East-West Direction"),
                 yaxis = list(title = "North-South Direction"),
                showlegend = T
        ) %>%
plotly::animation_opts(frame = 120,
                       transition = 5,
                       redraw = F)